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I.    Introduction 

This  report  describes  three  procedures  for  estimating  the 
position  of  an  object.   Each  procedure  is  based  on  a  model  that 
describes  the  position  of  an  object  in  terms  of  the  position  of  a 
point  on  a  plane  surface. 

The  first  procedure  is  developed  in  Section  II.   The 
procedure  is  a  method  for  combining  into  a  composite  estimate  two 
or  more  point  estimates  of  the  position  of  an  object  that  each 
have  confidence  regions  that  are  bounded  by  an  ellipse.   In  the 
model  on  which  it  is  based,  the  coordinates  of  each  estimate  are 
the  values  of  random  variables  that  have  a  bivariate  normal 
distribution  with  a  mean  vector  whose  elements  are  the  unknown 
coordinates  of  the  object.   A  user  of  the  procedure  must  specify 
the  coordinates  of  each  estimate  and  the  major  and  minor  axes  of 
the  ellipse  that  bounds  its  confidence  region. 

The  second  procedure  is  developed  in  Section  III.   The 
procedure  is  a  method  for  generating  a  position  estimate  from 
bearings  that  are  on  or  from  two  or  more  reference  points  of 
known  position.   In  the  model  on  which  it  is  based,  bearing 
errors  are  independent  normally  distributed  random  variables. 
Based  on  the  model,  the  coordinates  of  a  position  estimate  that 
is  determined  by  the  procedure  are  values  of  random  variables 
that  have  a  bivariate  normal  distribution  with  a  mean  vector 
whose  elements  are  the  unknown  coordinates  of  the  object.   A  user 
of  the  procedure  must  specify  the  standard  deviations  of  the 
bearing  errors.   An  extension  of  the  procedure  is  also  developed 


in  this  section.   It  is  a  method  for  generating  a  position 
estimate  from  bearings  that  are  on  or  from  two  or  more  reference 
points  for  which  the  position  of  one  or  more  of  the  reference 
points  is  uncertain. 

The  third  procedure  is  described  in  Section  IV.   The 
procedure  is  a  method  for  generating  a  position  estimate  that  is 
determined  from  lines  of  position.   In  the  model  on  which  the 
development  is  based,  lines  of  position  are  straight  lines,  lines 
of  position  that  are  determined  by  observation  are  parallel  to 
true  lines  of  position  and  the  algebraic  distance  between  a  line 
of  position  based  on  an  observation  and  a  true  line  of  position 
is  the  value  of  a  normally  distributed  random  variable.   A  user 
of  the  procedure  must  specify  the  standard  deviations  of  these 
random  variables. 

In  the  models  on  which  the  three  procedures  are  based, 
systematic  errors  represent  mean  values  that  are  set  to  zero. 
This  suggests  that  the  procedures  should  be  used  only  were  the 
bias  in  the  measurements  that  determine  bearings  or  lines  of 
position  is  known  or  is  known  to  be  negligible. 


II.   A  Composite  Position  Estimate 

The  procedure  that  is  developed  in  this  section  is  for 
combining  position  estimates  for  an  object  that  are  generated 
from  data  from  independent  sources.   It  is  based  on  the  following 
model:   In  a  rectangular  coordinate  system,  the  coordinates  of 
each  position  estimate  are  values  of  random  variables  that  are 
determined  by  an  independent  bivariate  normal  distribution  with  a 
covariance  matrix  whose  elements  are  known  and  with  a  mean  vector 
whose  elements  are  the  unknown  coordinates  of  the  object.    The 
above  model  implies  that  the  natural  logarithm  of  the  likelihood 
function   L   for  a  set  of   n   estimates  can  be  expressed  as 
follows:   log  L  =  K  -  §  2"  (xi  -  x)  2'1  (x^  -  x)   where   K   is  a 
constant,   xi  is  an  estimate  of  the  mean  vector  with  components 
xi  and  yi  ,   x  is  the  mean  vector  with  components   x   and   y,  the 
unknown  coordinates  of  the  object,  and   2,  is  the  covariance 
matrix  with  elements  o -  ,  o\-    and   0^.   The  composite  maximum 
likelihood  estimates   x   and  y   of  the  unknown  coordinates   x 
and  y  make   log  L  a  maximum  which  implies  the  following 
relations: 


d (log  L) 


=  0        and        d flog  L) 
x=x  dy 

y=y 


=  0 
x=x 

y=y 


After  evaluating  the  derivatives,  they  can  be  written  as 
A  x  +  B  y  =  D 
B  x  +  C  y  =  E 


where: 

A   =   Zj   aff       B  -  J}  b,    ,       C  -  2;   Cl  ,       D  =   2;    (a,xl    +   biYi)     , 

E   =   2;    (b1xi    +   Ciy.)     ,      a,   =      a2-   /   di  ,      b,   =   -   ffiX-   /   d1    , 

ci    =      a?i  /   <*i    >      di    "      a?x   *ij   "   aixy 
The   solutions   to   the   equations   are: 

x   =    (CD    -    BE)    /    (AC    -    B2) 

y   =    (AE    -    BD)    /    (AC    -    B2) 

Since  the  estimates   x   and   y   are  linear  combinations  of  the 
estimates   x1   and   yi  ,  they  are  determined  by  a  bivariate 
normal  distribution  that  is  defined  by  its  mean  vector  with 
components   x  and   y   and  its  covariance  matrix  with  elements 
a-,  a-   and  a-.      With   E   the  expected  value  operator: 

x  =  E{(CD  -  BE)/(AC  -  B2)  } 

y  =  E{(AE  -  BD)/(AC  -  B2)  } 

o\    =    E{(CD  -  BE)  -  E(CD  -  BE)}2/ (AC  -  B2)2 

=  (C2(F+I+2L)  -  2CB(H+K+M)  +  B2  (G+J+2N)  }/ (AC  -  B2)2 

ff|  =  E{(AE  -  BD)  -  E(AE  -  BD)}2/(AC  -  B2)2 

=    (A2(G+J+2N)     -    2AB(H+K+M)    +    B2  (F+I+2L)  }/ (AC    -    B2)2 

a-    =    E( [ (CD-BE) -E (CD-BE) ][ (AE-BD) -E (AE-BD) ] )/ (AC    -    B2)2 


xy 


=    (  (AC+B2)  (H+K+M)-CB(F+I+2L)-BA(G+J+2N)  }/(AC    -    B2)2 


where: 

F  =  2?   a2a2-    ,      G  =  2?  b2a2-    ,      H  =   2?   aba2-    .      1  =  2?  b2a2-    , 

J  =  SJ  c2a2y    ,      K  =  2J  blC,a2y    ,      L  =  Zj   a^a,--    ,      M  =  2J    (a^+b2)  a,-y 

Without  loss  of  generality,  the  origin  of  the  coordinate 
system  can  be  the  point  determined  by  the  composite  estimate,  its 
positive  y-axis  can  be  north  and  its  positive   x-axis  east.   In 
Section  II   of  this  report,  it  is  shown  that  the  axes  of  the 
ellipse  that  bounds  a  minimum  area  confidence  region  for  a 
composite  estimate  are  coincident  with  a  primed  coordinate  system 
which  is  transformed  to  this  coordinate  system  by  a  coordinate 
axes  rotation  equal  to   7   where   7   is  defined  by  the  relation: 
tan  27  =  2ct;-  /(a?  -  a\)    . 

It  can  be  shown  that  a  minimum  area  confidence  region  for  a 
two  dimensional  estimator  whose  components  have  a  bivariate 
normal  distribution  is  bounded  by  an  ellipse  whose  major  axes  are 
determined  by   2ka-,   and   2ka~.   where  k  =  [-2  log(l  -  p)]1/2 

x  y 

and   x'  and   y1   refer  to  a  coordinate  system  whose  axes  are 
parallel  to  the  axes  of  the  bounding  ellipse.   In  this  system: 

a\.   =     o\   cos2  7-2  o-   sin  7  cos  7   +   a?  sin2  7 

a2.  =  a\   sin2  7   +   2  o-~   sin  7  cos  7   +   o\   cos2  7 

y  x  xy  y 

a-.-.   =      0 

x  y 


The  model  that  is  described  above  can  be  used  to  find  a 
composite  position  estimate  from  two  or  more  position  estimates 
for  an  object.   To  show  this,  first  assume  that  the  estimates 
satisfy  the  conditions  that  are  required  by  the  model.   Then,  let 
6i   be  the  direction  of  the  major  axis  of  the  ellipse  that  bounds 
the  confidence  region  that  is  associated  with  the  i  th   position 
estimate.   Next,  orient  the  primed  coordinate  system  defined 
above  so  that  its   x'-axis  is  coincident  with  the  major  axis  of 
the  ellipse.   Then,  the  elements  of  the  covariance  matrix  of  the 
bivariate  normal  distribution  that  determines  the  estimate  are 
given  by  the  following  relations: 

az~   =  or2.,  sin2  6       +      a2-.  cos2  8 

IX        IX   ^j-*1    «*  i         v   ly   WB    "l 

a2-  =  a2-,  cos2  6i  +  a2~.  sin2  5i 

aiiy  =  (aix-   "  aV>  sin  6i  cos  5i 

As  an  example  of  the  use  of  the  procedure,  suppose  that  each 
bounding  ellipse  were  a  circle.   In  this  case,  the  coordinates  of 
each  position  estimate  would  be  determined  by  a  circular  normal 
distribution  and,  for  i  =  1  to  n  ,  the  elements  of  the  covariance 
matrices  of  the  circular  normal  distributions  would  satisfy  the 
relations:   ai;  =  a-  =  ai   and  a--  =  0.   The  composite  estimate 
would  be: 

x  =  [2J  (i,  /  o,)]   /   [2;  (l  /  a,)] 

y  =  [2;  (y]  /   a,)]  /  [  (2J  l  /  oj  ] 


And  the  elements  of  the  covariance  matrix  of  the  distribution 
that  determine  the  coordinates  of  the  estimate  would  be: 

a\   -  n  /  [2;  (1  /  o,)f    ,  a]   =  n  /[(Sj  (1  /  a,)]2   and   ax-  =  0  . 

Consequently,  the  distribution  would  also  be  circular  normal  and 
the  bounding  curves  for  minimum  area  confidence  regions  would  be 
circles. 

A  program  called  COMP  that  is  based  on  the  procedure  that  is 
developed  in  this  section  is  listed  in  Appendix  1.   COMP  can  be 
used  to  generate  composite  position  estimates  and  corresponding 
confidence  regions.   As  an  example  of  its  use,  suppose  three 
independent  position  estimates  are  given  in  Table  1. 

TABLE  1 

i      xi        yi      6i       MJ1        MIi        k 

1  -3.7       18.1      59°       36  20         2 

2  11.8         8.4     105°        38  10  2 

3  0  0      146°        50  24  2 

In  Table  1,   MJ   is  the  length  of  the  major  axis  and  MI   is  the 
length  of  the  minor  axis  of  the  ellipse  bounding  a  minimum  area 
confidence  region.   If  the  distance  unit  is  the  nautical  mile, 
the  coordinates  of  the  composite  estimates  of  the  coordinates  of 
the  object  are:   x  =  -2.69  nm   and  y  =  12.41  nm.   For   k  =  2, 
the  lengths  of  the  major  and  minor  axes  of  the  bounding  ellipse 
of  the  minimum  area  confidence  region  are:   MJ  =  17.3  3  nm  and 
MI  =  8.85  nm  and  the  direction  of  its  major  axis  is   103.77°. 
The  results  are  shown  in  Figure  1. 
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Figure  1.   The  position  estimates  and  their  associated  confidence 
regions  for  the  example  on  Page  7.   The  position  estimates  are  at 
the  center  of  their  corresponding  confidence  region.   The  numbers 
indicate  the  order  of  the  estimates  in  the  table  on  Page  7.   The 
composite  confidence  region  is  in  black. 


III.  A  Position  Estimate  Based  on  Bearings 

The  procedure  that  is  developed  in  this  section  is  for 
determining  a  position  estimate  for  an  object  from  a  set  of 
observed  bearings  from  or  on  two  or  more  reference  points  of 
known  positions.   In  the  development,  reference  points  are 
referred  to  as  stations.   The  procedure  is  based  on  the  following 
model:   A  set  of  bearings  taken  on  an  object  from  a  station  are 
the  values  of  independent  normally  distributed  random  variables 
with  known  standard  deviations  but  with  means  that  are  equal  to 
the  unknown  true  bearings.   The  model  implies  that  the  natural 
logarithm  of  the  likelihood  function   L   for  a  set  of   n 
bearings  taken  on  an  object  can  be  expressed  as  follows: 
log  L  =  K  -  1/2  2J  (^  -  <f>y)2  /   &]        where   K   is  a  constant,  0, 
is  an  observed  bearing,  0i   is  the  unknown  true  bearing  and   ei 
is  the  known  standard  deviation.   For   i  =  1  to  n  ,  the  maximum 
likelihood  estimates  <f>.      are  the  solutions  of  the   n   equations: 


a (loo  L) 


=       0 


subject  to  the  constraint  that  the  bearing  lines  determined  by 
the  estimates  must  all  pass  through  a  common  point.   Because  of 
this  constraint,  the  number  of  independent  estimates  is  reduced 
from  n   to   2.   This  is  equivalent  to  stating  that  any  two  of 
the  estimates  determine  a  common  point.   The  common  point  is  the 
maximum  likelihood  estimate  for  the  position  of  the  object. 


In  order  to  impose  the  constraint  on  the  estimates,  the 
likelihood  function  will  be  reformulated  in  a  rectangular 
coordinate  system.   To  do  this,  first,  consider  the  following 
definitions: 

u,  =  r,  (^  -  0,)  ,   vi  =  ri  (0,  -  fis)      and  w,  =  r,  (*,  -  0,) 

where   i  =  1  to  n   is  the  observation  number  and  the  observations 
are  from  two  or  more  stations.   In  these  definitions,   j3i   is  the 
bearing  of  an  initial  estimate  of  the  position  of  an  object  and 
ri   is  its  range  from  the  station  associated  with  the   i  th 
observation.   It  is  also  the  radius  of  a  circle  that  is  centered 
on  the  station  and  passes  through  the  initial  estimate  as  shown 
in  Figure  2.   With  0S    ,  <f>-      and   /^   in  radians,  ui  ,  vi   and  wi 
are  the  length  of  arcs  on  this  circle.   And,   ui  =  wi  -  vi  .   In 
this  relation,  vi   involves  the  unknown  true  bearing   0i  ,  but  wi 
is  known. 

In  order  to  express   vi   in  a  convenient  form,  consider  a 
rectangular  coordinate  system  whose  origin  is  at  the  initial 
estimate,  whose  positive   y-axis  is  north  and  whose  positive 
x-axis  is  east.   With   x   and   y   the  coordinates  of  the  object, 
to  first  order:   vi  =  x  cos  ($■      -  y  sin  y3i  .   This  implies: 

0i  =  /3i  +   (x  cos  /3i  -  y  sin  0J    /   ri 

And,  to  first  order,  this  expresses  the  constraint  on  the   0i 

and  on  their  maximum  likelihood  estimates  and  it  implies  that  the 

bearing  lines  determined  by  /^  ,  0i   and  8-      are  parallel. 
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Figure  2.   In  the  program  that  implements  the  procedure,  the 
origin  of  the  xy-coordinate  system  is  coincident  with  the  common 
reference  point  rather  than  with  the  initial  estimate  and  the 
coordinates  of  the  initial  estimate  are   x*   and   y*.   The 
program  is  discussed  later  in  the  section. 
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In  order  to  determine  the  estimates   0i   to  first  order,  use 
the  above  relation  to  replace   0i   by  the  known  values  fl-      and 
the  two  unknown  values   x   and  y.   Next,  determine  the  maximum 
likelihood  estimates   x   and   y   for  the  unknown  coordinates   x 
and   y  which,  to  first  order,  are  the  unknown  coordinates  of  the 
object.   That  is,  find  the  solutions  to  the  equations  determined 
by: 


d flog  L) 

dx 


=  0 


and 


x=x 

y=y 


a  flog  L) 
dy 


=  0 


x=x 

y=y 


Evaluating  the  derivatives,  the  equations  are: 

Zj    [  (wi      -      x   cos   fis      +      y   sin   $s)    cos   /Jj    /    a,   =    0 
2"    [  (wi      -      x  cos   ^i      +      y  sin  fis)    sin  fl  s]    /   ai   =   0 

where      ai    =   ri    eia      These   equations   can   also   be  written   as 

A   x   +   B  y   =   D 
B   x  +   C  y   =   E 

and  their  solutions  as: 

x   =    (CD   -   BE)    /    (AC   -    B2) 
y   =    (AE    -    BD)    /     (AC    -    B2) 

where 


A 

C 
E 


2;    (cos2  fij    /   o\    ,      B  =  2;    (cos   fis   sin  fis)    /   o\    , 
SJ    (sin2  0,)   /  a]   ,      D  =  Sj    (Wi   cos  fis)    /  o\ 
2j    (wi    sin   fit)    /   o\    . 
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If  the  distribution  of  the  random  variables  that  determine   x 
and   y   have  a  known  distribution,  then,  in  principle,  a 
confidence  region  for  an  estimate  can  be  constructed.   With  this 
in  mind,  the  distribution  of  the  estimators  that  determine   x 
and   y   will  be  determined  next.   If   X   and   Y   represent  the 
estimators,  that  is,  the  random  variables  that  determine   x   and 
y  ,  then, 

X  =  2j  (Wi  /  a2)     (B  sin  ^   -   C  cos  fi .)    /    (B2  -  AC) 

Y  =  2j  (Wi  /  O2)     (A  sin  0.      -   B  cos  fi .)    /    (B2  -  AC) 

where  Wi  =  ri  (0i  -  fl })     .   Since  the   8i   are  normally  distributed 
random  variables,   the  Wi   are  also.   And,  since  the  estimators   X 
and   Y   are  linear  combinations  of  the   Wi  ,   they  have  a 
bivariate  normal  distribution. 

If  fl s   =  0i  ,  for  i  =  1  to  n,  that  is,  if  the  initial 
estimate  is  the  position  of  an  object  then,  with   E   again  the 
expected  value  operator,  E(W^)  =  0  ,  since  Wi  =  ri  (0i  -  fi ,)     .  In 
this  case: 

o\  =   2J  (1/a2,)  (B  sin  j81  -  C  cos  ^)2/(B2  -  AC)2 

a]   =  2J  (l/a2i)(A  sin  ^   -   B  cos  P,)2/(B2   -   AC)2 

a--  =  2J  (1/ff2,)  (B  sin  ^  -  C  cos  ^,)  (A  sin  >3  i  -  B  cos  ^)  /  (B2   -   AC)2 

Using  the  definitions  for  A  ,  B  ,  and   C  ,   that  are  given 
above,  these  equations  can  be  written  as: 
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a\   =  C  /  (AC  -  B2) 

o]    =    A    /     (AC  -  B2) 

a--  =  B  /  (AC  -  B2) 

To  the  degree  that  the  initial  estimate  is  near  the  position  of 
an  object,  these  equations  define  approximations  for  the  elements 
of  the  covariance  matrix  of  the  estimators   X   and   Y. 

Since,  in  general,  a—   will  not  be  zero,  the  axes  of  the 
ellipse  that  bounds  a  minimum  area  confidence  region  for   x   and 
y  will  not  be  parallel  to  the   xy-coordinate  axes.   To  determine 
the  orientation,  translate  the   xy-coordinate  system  so  that  its 
origin  is  at  the  point  determined  by  x   and   y  .   Then,  rotate 
this  coordinate  system  until  its  axes  coincide  with  the  axes  of 
the  bounding  ellipse.   Next,  represent  coordinates  in  this  system 
by   x'   and   y'.   If  the  rotation  angle  is   7  ,  then  the 
coordinates  of  a  point  in  the  two  systems  are  related  by  the 
following  general  transformation  equations: 

x'  =  x  cos  7   -  y  sin  7 
y'  =  x  sin  7   -   y  cos  7 

These  transformation  equations  can  be  used  to  determine  the 
elements  of  the  covariance  matrix  of  the  distribution  in  the 
primed  coordinate  system  from  the  covariance  elements  of  the 
distribution  matrix  in  an  unprimed  coordinate  system,  by  using 
the  following  equations: 
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a-.  =  ok   cos  7   -   2  a--  sin  7  cos  7   +    a?  sin  7 
a?.  =   a?  sin2  7   +   2  a;-  sin  7  cos  7   +   a?  cos2  7 

ax'y'  =   (ai  "  °\)    sin  7  cos  7   +  oky    (cos2  7  -  sin2  7) 

Since   X'  and   Y1  are  independent  random  variables  when  the  axes 
of  the  coordinate  system  are  parallel  to  the  axes  of  the  bounding 
ellipse,   <7j.-.  =  0   which  implies  that   7  satisfies  the  equation: 

tan  27  =  2o~   /{o\   -      o\)     . 

If  the  position  of  an  object  were  known,  and  the  estimated 
position  were  at  the  true  position,  then  the  probability   p   that 
an  estimated  position  would  be  in  a  region  bounded  by  an  ellipse 
with  axes   2ka-.   and   2ka-.  that  is  centered  on  the  true  position 
is   given  by   p  =  1  -  exp(-k2/2).   This  result  can  be  obtained  by 
integrating  over  the  region  bounded  by  the  ellipse  the  bivariate 
normal  density  that  determines   X'  and   Y1.   Note  that   k   is  a 
measure  of  the  size  of  the  region  bounded  by  the  ellipse  and,  for 
a  probability  of  containment   p  ,  the  corresponding  value  of   k 
is  determined  by   k  =  [-2  log(l  -  p))]1/2  •    These  relations 
provide  a  basis  for  determining  confidence  regions  for  the 
estimates   x   and   y   of  the  position  of  an  object  that  are  given 
above . 

A  program  called  PEST  that  is  listed  in  Appendix  2  can  be 
used  to  generate  position  estimates  and  corresponding  confidence 
regions  from  bearing  observations  on  or  from  two  or  more 
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stations.   For  a  bearing  on  a  station,  9      is  the  reciprocal  of 
the  observed  bearing.   In  the  program,  the  initial  estimate  is 
the  point  of  intersection  of  the  two  bearing  lines  that  are 
determined  by  the  first  two  entered  bearings  with  each  from  a 
different  station.   After  the  first  two  bearings  have  been 
entered,  the  program  determines  the  coordinates   x*   and  y*   of 
the  initial  estimate  by  solving  the  following  equations: 

x*  sin  (92   -   9:)    =  [pl   sin  (c^  -  9  {)  ]  sin  92 

-  [p2   sin  (a2   -  62)  ]  sin  81 

y*    sin    (92   -    9 j)    =    [p1    sin    (a1    -    8  x)  ]    cos    62 

-  [p2   sin    (a2   -   9 2)  ]    cos   9X 

where  al        and  pl        are   the  bearing  and  range  of  the  first 
station  and  a2      and  p2      are  the  bearing  and  range  of  the  sscond 
station  from  a  common  reference  point.   This  method  of  finding 
the  initial  estimate  is  based  on  one  that  is  used  in  Reference  3. 
The  point  estimates  a      and  p      of  the  bearing  and  range  of  an 
object  from  a  common  reference  point  are  shown  in  Figure  3. 
As  an  example,  consider  the  data  in  Table  2. 


TABLE  2 

Station 

Observed 

Station 

Station 

Bearing 

Number 

Bearing 

Bearing 

Range 

Error 

1 

38° 

334° 

13500 

4° 

2 

324° 

50° 

11350 

3° 

3 

3° 

0° 

0 

4° 
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East 


common 
reference  point 


Figure  3.   A  confidence  region  bounded  by  an  ellipse  and  a  primed 


coordinate  system  in  which 


x  y 


is  zero 


The  center  of  the 


ellipse  and  the  origin  of  the   primed  and  the  unprimed  coordinate 

systems  are  at  the  estimated  position  of  an  object.   The 

direction  of  the  major  axis  of  the  ellipse  is  8      where 

0°  <  6  <  180°.   The  rotation  angle  from  the  unprimed  to  the 

primed  coordinate  system  is   7  .   The  estimates  a      and  p      are 

the  bearing  and  range  of  the  object  from  a  common  reference 

point. 
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The  order  of  data  entry  in  the  program  is  determined  by  the 
station  number.   Note  that  the  reference  point  is  at  the  location 
of  Station  3.   If  the  distance  is  in  meters,  then,  for  the  three 
observations,  the  estimates  generated  by  the  program  for  the 
bearing  and  range  of  the  object  from  Station  3  are:  a   =  0.15° 
and  p   =  19554  meters.   For   k  =  2,  the  lengths  of  the  major  and 
minor  axes  of  the  ellipse  bounding  the  minimum  area  confidence 
region  are:   MJ  =  3426  meters   and   MI  =  2260  meters   and  the 
direction  of  its  major  axis  is   12.48°. 

An  extension  of  the  above  procedure  is  developed  next.   The 
extension  accounts  for  cases  where  not  all  of  the  positions  of 
the  reference  stations  can  be  considered  to  be  known.   In  these 
cases,  if  the  station  position  were  to  become  known,  then,  for  a 
bearing  line  that  is  determined  by  an  observation,  the  bearing 
line  that  is  drawn  from  an  assumed  station  position  will  be 
parallel  to  the  bearing  line  that  is  drawn  from  the  station 
position.   The  procedure  is  based  on  this  relation  and  on  the 
model  with  the  addition  of  the  following  condition:   The 
coordinates  of  an  uncertain  station  position  are  determined  by  a 
distribution  which  is  bivariate  normal  with  a  covariance  matrix 
whose  elements  are  known  and  a  mean  vector  whose  elements  are  the 
unknown  coordinates  of  the  station.   This  implies  that  for  a 
bearing  determined  by  an  observation,  the  algebraic  distance 
between  the  bearing  line  from  an  assumed  station  position  and  the 
bearing  line  from  the  station  position  is  a  random  variable   S 
that  has  a  normal  distribution  with  mean  zero  and  a  standard 
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deviation   as   that  is  determined  by  the  elements  of  the  known 
covariance  matrix. 

The  geometry  of  the  model  is  shown  in  Figure  4.   There,  8 
is  an  observed  bearing  and  the  common  origin  of  the  coordinate 
systems  is  the  assumed  station  position.   For  the  unprimed 
system,  the  positive   y-axis   direction  is  north  and  the  positive 
x-axis   direction  is  east.   For  the  primed  system,  the   x'-axis 
is  coincident  with  the  major  axes  of  an  ellipse  that  bounds  a 
confidence  region  for  the  assumed  position.   The  direction  of  the 
major  axis  is   6   where   0°  <  6  <  180°.    For  the  double  primed 
system,  the   y"-axis  is  coincident  with  the  assumed  bearing  line. 

The  effect  of  station  position  uncertainty  is  accounted  for 
by  the  algebraic  distance   s   between  a  bearing  line  from  an 
assumed  station  position  and  the  parallel  bearing  line  from  the 
station  position.   In  effect,  to  first  order,  the  following 
bearing  lines  are  also  parallel:  the  bearing  line  from  a  station 
to  the  initial  estimate,  the  bearing  line  from  the  station  to  an 
object  and  a  bearing  line  determined  by  an  observation.   Hence, 
the  random  variable   Ui  =  ri  (8,  -  <j> f)   is  the  algebraic  distance 
between  the  observed  bearing  line  and  the  bearing  line  from  the 
station  to  an  object.   And  Uei  =  ri  (8,  -  0,)  ±  S   is  a  random 
variable  that  determines  the  algebraic  distance  between  the 
observed  bearing  line  from  an  uncertain  station  position  and  the 
bearing  line  from  the  station  position  to  the  object.   Since  U, 
and   S   are  independent  random  variables,  the  standard  deviation 
of  Uei  is:   aei  =  [a?  +  o\]m    . 
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In  the  extended  procedure,  for  a  station  of  uncertain 
position,  the  standard  deviation   ai  =  r1ei   is  replaced  by: 

o   .  =  [o'< -   +   a2]172 

ei     !■  i      s  J 

where 

o\   =   a2.  sin2  (6-6)      +   a2.  cos2  (6  -  6) 

based  on  Figure  4. 

In  Figure  4,  note  that   kcr2.   and   ka2.   are  the  lengths  of 

x  y 

the  major  and  minor  axes  of  the  bounding  ellipse  for  some  value 
of  k.   For  an  uncertain  station  position,   ^i   the  bearing  and 
ri   the  range  of  the  initial  estimate  from  the  station  position 
are  determined  by  the  assumed  station  position.   These 
relationships  are  illustrated  in  Figure  5. 

In  the  development,  although  the  position  of  one  or  more  of 
the  stations  may  be  unknown,  the  position  of  a  common  reference 
point  is  known.   This  implies  that  the  assumed  station  positions 
and  the  position  of  the  initial  estimate  are  also  known. 
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Figure  5.   In  the  development  for  the  extended  procedure,  fi      and 
r   are  the  range  and  bearing  of  the  initial  estimate  from  the 
assumed  station  position. 


22 


IV.   A  Position  Estimate  Based  on  Lines  of  Position 

The  procedure  that  is  described  in  this  section  is  for 
determining  a  position  estimate  for  an  object  from  two  or  more 
lines  of  position  that  can  be  considered  to  be  straight  lines  in 
the  neighborhood  of  the  object.   The  procedure  is  based  on  the 
following  model:   A  line  of  position  is  a  straight  line.   A  line 
of  position  that  is  based  on  an  observation  is  parallel  to  the 
true  line  of  position.   The  algebraic  distance   u   between  a  line 
of  position  that  is  based  on  an  observation  and  the  true  line  of 
position  is  the  value  of  a  normally  distributed  random  variable 
with  zero  mean  and  known  standard  deviation. 

Labeling  lines  of  position  by  the  index   i  ,  the  model 
implies  that  the  algebraic  distance   wi   between  the  i  th  line  of 
position  based  on  an  observation  and  a  parallel  line  through  an 
initial  estimate  of  an  unknown  position  is  related  to  the 
algebraic  distance   ui   between  it  and  the  true  line  of  position 
by:   ui  =  w^-  vi   where  vi   is  the  algebraic  distance  between  the 
i  th  true  line  of  position  and  the  parallel  line  through  the 
initial  estimate.   Now,  consider  a  rectangular  coordinate  system 
whose  origin  is  at  the  initial  estimate,  whose  positive   y-axis 
is  north  and  whose  positive   x-axis  is  east.   With   x   and   y 
the  coordinates  of  the  unknown  position,  to  first  order: 

vi  -   x  cos  ^i   -  y  sin  /^ 

where  fls      is  the  direction  of  the  i  th  true  line  of  position  and 
0°  <  ($■    <    180°.   By  referring  to  the  development  in  Section  III, 

23 


note  that  the  estimates   x   and   y  that  are  developed  there  also 
apply  here.   However,  in  this  development  the  standard  deviation 
o   of  the  random  variable  that  determines  the  distance   u  must  be 
known  for  each  of  the  observed  bearing  lines. 

Figure  6   illustrates  the  geometry  associated  with  the 
development. 
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Figure  6.   The  geometry  for  three  lines  of  position 


The 


direction  of  the  lines  are   /} 


f}?      and   /L.   With  the  initial 


estimate  of  the  unknown  position  at  the  intersection  of  line  1 
and  line  2  ,  both  Wj  and  w2  are  equal  to  zero.  This  implies 
that  the  final  estimate  of  the  position  and  the  initial  estimate 
of  the  position  are  equal.  If  the  three  lines  of  position  were 
obtained  by  sextant  observations,  the  assumed  position  should  be 
the  reference  position.  In  this  case,  the  angles  al  ,  a~  and 
a3   would  be  the  azimuth  angles. 
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Appendix  1.  The  Program  COMP 

COMP  is  a  program  with  which  to  implement  the  procedure 
developed  in  Section  II   for  generating  a  composite  estimate  of 
the  position  of  an  object  from  two  or  more  independent  estimates 
of  the  position  of  the  object.   The  independent  estimates  must  be 
specified  in  terms  of  coordinates  in  a  rectangular  coordinate 
system  whose  positive   y-axis  is  north  and  positive   x-axis  is 
east.   In  addition,  a  confidence  region  for  each  estimate  must  be 
specified  in  terms  of  the  direction   6   of  the  major  axis  of  a 
bounding  ellipse  and  either  its  size   k  or  its  containment 
probability   p   (confidence) . 

The  program  listing  follows: 

10  CLS  :  PI  =  4  *  ATN(l) 

20  INPUT  "print  data  yes/no  (y/n)";  A$ 

30  IF  A$  =  "N"  OR  A$  =  "n"  THEN  P$  =  "N":  GOTO  60 

40  IF  A$  =  "Y"  OR  A$  =  "y"  THEN  P$  =  "Y"  ELSE  GOTO  20 

50  LPRINT  :  LPRINT  :  LPRINT  :  LPRINT  "Composite  Position  Estimation  Program:  COMP.BAS" 

60  PRINT  :  INPUT  "estimates  to  be  combined";  NO 

70  IF  P$  =  "N"  THEN  GOTO  90 

80  LPRINT  :  LPRINT  "estimates  to  be  combined    ",  NO 

90  INPUT  "elliptical  confidence  regions  by  size  or  probability  (s/p)";  F$ 

100  IF  F$  =  "S"  OR  F$  =  "s"  OR  F$  =  "P"  OR  F$  =  "p"  THEN  GOTO  110  ELSE  GOTO  90 

110  IF  F$  =  "P"  OR  F$  =  "p"  THEN  B$  =  "probability"  ELSE  B$  =  "size" 

120  FOR  10  =  1  TO  NO 

130  PRINT  :  INPUT  "position  estimate  x-coordinate";  X 

140  INPUT  "position  estimate  y-coordinate";  Y 

150  INPUT  "containment  ellipse  major  axis  direction";  DO 

160  IF  DO  >=  180  OR  DO  <  0  THEN  GOTO  150 

170  INPUT  "ellipse  major  axis  length";  MA:  SMA  =  MA  /  2 

180  INPUT  "ellipse  minor  axis  length";  MI:  SMI  =  MI  /  2 

190  IF  F$  =  "P"  OR  F$  =  "p"  THEN  GOTO  220 

200  INPUT  "containment  ellipse  size";  SO:  P0  =  1  -  EXP(-S0  *  SO  /  2) 

210  PRINT  "containment  probability  =  ";  P0:  GOTO  240 

220  INPUT  "containment  probability";  P0:  SO  =  SQR(-2  *  L0G(1  -  P0)) 

230  PRINT  "containment  ellipse  size  =  ";  SO 

240  IF  P$  =  "N"  THEN  GOTO  320 

250  LPRINT  :  LPRINT  :  LPRINT  "position  estimate  x-coordinate",  X 

260  LPRINT  "position  estimate  y-coordinate",  Y 

270  LPRINT  "ellipse  major  axis  direction",  DO 

280  LPRINT  "ellipse  major  axis  length   ",  MA 

290  LPRINT  "ellipse  minor  axis  length   ",  MI 

300  LPRINT  "ellipse  size  ",  SO 

310  LPRINT  "containment  probability     ",  P0 

320  SX  =  SMA  /  SO:  SY  =  SMI  /  SO 

330  DO  =  DO  *  PI  /  180 

340  U  =  SX  *  SX  *  SIN(DO)  *  SIN(DO)  +  SY  *  SY  *  C0S(D0)  *  C0S(D0) 

350  V  =  SX  *  SX  *  C0S(D0)  *  COS(DO)  +  SY  *  SY  *  SIN(DO)  *  SIN(DO) 

360  W  =  (SX  *  SX  -  SY  *  SY)  *  SIN(DO)  *  C0S(D0) 

370Z=U*V-W*W:A=A+V/Z:B=B-W/Z:C=C+U/Z:D=D+V/Z*X-W/Z*Y 
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380E=E-W/Z*X+U/Z*Y:F=F+V/Z*V/Z*U:G=G+W/Z*W/Z*U 

390H=H-V/Z*W/Z*U:I=I+W/Z*W/Z*V:J=J+U/Z*U/Z*V 

400K=K-W/Z*U/Z*V:    L   =   L-V/Z*W/Z*W 

410M  =  M+(V/Z*U/Z  +  W/Z*W/Z)*W:N  =   N-W/Z*U/Z*W 

420  NEXT  10 

430  Q  =  A  *  C  -  B  *  B:  X  =  (C  *  D  -  B  *  E)  /  Q:  Y  =  (A  *  E  -  B  *  D)  /  Q 

440R=(C*C*(F+I+2*L)-2*C*B*(H+K+M)+B*B*(G+J+2*N))/(Q*Q) 

450S=(A*A*(G+J+2*N)-2*A*B*(H+K+M)+B*B*(F+I+2*l))/(Q*Q) 

460  T  =  ((A  *  C  +  B  *  B)  *  (H  +  K   +  M)  -  C  *  B  *  (F  +  I  +  2  *  L)  -  B  *  A  *  (G  +  J  +  2  *  N))  /  (Q  *  Q) 

470  CLS  :  PRINT  "composite  estimate  x-coordinate  =  ";  X 

480  PRINT  "composite  estimate  y-coordinate  =  ";  Y 

490  IF  P$  =  "N"  THEN  GOTO  520 

500  LPRINT  :  LPRINT  :  LPRINT  :  LPRINT  "composite  estimate  x-coordinate",  X 

510  LPRINT  "composite  estimate  y-coordinate",  Y 

520  IF  S  =  R  THEN  A  =  0:  GOTO  600:  REM  a  circle 

530  A  =  .5  *  ATN(2  *  T  /  (S  -  R)) 

540  Rl  =  R  *  COS(A)  *  C0S(A)  -  2  *  T  *  C0S(A)  *  SIN(A)  +  S  *  SIN(A)  *  SIN(A) 

550  SI  =  R  *  SIN(A)  *  SIN(A)  +  2  *  T  *  COS(A)  *  SIN(A)  +  S  *  COS(A)  *  COS(A) 

560  A  =  A  *  180  /  PI 

570  IF  Rl  >  SI  THEN  A  =  A  +  90:  GOTO  590 

580  RT  =  Rl:  Rl  =  SI:  SI  =  RT 

590  IF  A  <  0  THEN  A  =  A  +  180 

600  PRINT  :  PRINT  :  INPUT  "confidence  region  by  size  or  probability  (s/p)";  G$ 

610  IF  G$  =  "S"  OR  G$  =  "s"  THEN  GOTO  630 

620  IF  G$  =  "P"  OR  G$  =  "p"  THEN  GOTO  650  ELSE  GOTO  600 

630  PRINT  :  INPUT  "containment  ellipse  size";  SO:  PO  =  1  -  EXP(-S0  *  SO  /  2) 

640  PRINT  :  PRINT  "containment  probability  =  ";  PO:  GOTO  670 

650  PRINT  :  INPUT  "containment  probability";  PO:  SO  =  SQR(-2  *  LOG(l  -  PO)) 

660  PRINT  :  PRINT  "containment  ellipse  size 

670  PRINT  "ellipse  major  axis  length  =  ";  2 

680  PRINT  "ellipse  major  axis  direction  =  "; 

690  PRINT  "ellipse  minor  axis  length  =  ";  2 

700  IF  P$  =  "N"  THEN  GOTO  760 

710  LPRINT  :  LPRINT  "containment  ellipse  size 

720  LPRINT  "containment  probability 

730  LPRINT  "ellipse  major  axis  length 

740  LPRINT  "ellipse  major  axis  direction 

750  LPRINT  "ellipse  minor  axis  length 

760  PRINT  :  PRINT  :  INPUT  "continue  or  quit  (c/q)";  G$ 

770  IF  G$  =  "C"  OR  G$  =  "c"  THEN  GOTO  600 

780  IF  G$  =  "Q"  OR  G$  =  "q"  THEN  GOTO  790  ELSE  GOTO  770 

790  END 


SO 
A 
SO 

SO 

*  SQR(Rl) 

*  SQR(Sl) 

PO 
2  * 
A 

2  ' 

",  SO 

'   SO  *  SQR(R1 
'   SO  *  SQR(S1 
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Appendix  2.   The  Program  PEST 

PEST  is  a  program  with  which  to  implement  the  procedure 
developed  in  Section  III   for  determining  a  position  estimate  for 
an  object  from  a  set  of  observed  bearings  from  or  on  two  or  more 
reference  points  of  known  positions.   In  the  development,  the 
reference  points  are  referred  to  as  stations.   The  bearing  error 
associated  with  an  observation  must  be  specified  for  each 
station.   In  addition,  the  position  of  each  station  must  be 
specified  in  terms  of  its  bearing  and  range  from  a  single 
reference  point  which  may  be  one  of  the  reference  points.   The 
program  generates  an  estimate  for  the  position  of  an  object  and  a 
confidence  region  for  the  estimate  which  is  specified  in  terms  of 
the  direction  8      of  the  major  axis  of  a  bounding  ellipse  and 
either  its  size  k  or  its  containment  probability  p 
(confidence) . 

The  program  listing  follows: 

10  CLS  :  PI  =  4  *  ATN(l):  DIM  A(7):  I  =  0 

20  INPUT  "print  data  (y/n)";  A$ 

30  IF  A$  =  "N"  OR  A$  =  "n"  THEN  P$  =  "N":  GOTO  60 

40  IF  A$  =  "Y"  OR  A$  =  "y"  THEN  P$  =  "Y"  ELSE  GOTO  20 

50  LPRINT:  LPRINT:  LPRINT  :  LPRINT  "Position  Estimation  Program:  PEST.BAS" 

60  INPUT  "bearings  on  object  or  bearings  on  stations  (o/s)";  B$ 

70  IF  B$  =  "0"  OR  B$  =  "o"  THEN  K  =  0:  GOTO  120 

80  IF  B$  =  "S"  OR  B$  =  "s"  THEN  K  =  1  ELSE  GOTO  60 

90  IF  P$  =  "N"  THEN  GOTO  120 

100  IF  K  =  0  THEN  B$  =  "  object"  ELSE  B$  =  "  stations" 

110  LPRINT  :  LPRINT  "bearings  on"  +  B$:  LPRINT 

120  CLS  :  PRINT  "input  number  "  +  STR$(I  +  1) 

130  INPUT  "observed  bearing  in  decimal  degrees";  P 

140  INPUT  "station  bearing  in  decimal  degrees";  Q 

150  INPUT  "station  range";  R 

160  INPUT  "bearing  error  in  decimal  degrees";  BE 

170  IF  P$  =  "N"  THEN  GOTO  230 

180  LPRINT  :  LPRINT  "input  number  "  +  STR$(I  +  1) 

190  LPRINT  "observed  bearing  in  decimal  degrees",  P 

200  LPRINT  "station  bearing  in  decimal  degrees",  Q 

210  LPRINT  "station  range  in  chosen  units",  R 

220  LPRINT  "bearing  error  in  decimal  degrees",  BE 

230  IF  K  =  1  THEN  P  =  P  +  180 

240  P  =  P  *  PI  /  180:  Q  =  Q  *  PI  /  180:  BE  =  BE  *  PI  /  180 

250  IF  I  =  2  THEN  GOTO  350 

260  I  =  I  +  1:  A(I  -  1)  =  P:  A(I  +  1)  =  Q:  A(I  +  3)  =  R:  A(I  +  5)  =  BE 
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270 

IF  I  =  1 

THEN  GOTO  120 

280 

X  =  A(4) 

*  SIN(A(2)  -  A(0)):  Y  =  A(5)  *  SIN(A(3)  - 

*(!)) 

Z  =  SIN( 

A(0)) 

290 

IF  Z  =  0 

THEN  GOTO  900 

300 

U  =  (X  * 

SIN(A(1))  -  Y  *  SIN(A(0)))  /  Z:  V  =  (X  *  C0S(A(1))  -  Y 

*  COS(A(0)) 

/  z 

310 

FOR  J  = 

)  TO  1 

320 

P  =  A(J) 

Q  =  A(J  +  2):  R  =  A(J  +  4):  BE  =  A(J  +  6) 

:  GOSUB 

330 

NEXT  J 

340 

GOTO  360 

350 

GOSUB  940 

360 

CLS  :  INPUT  "continue  data  input  or  generate  an  est 

imate 

(c/g)"; 

C$ 

370 

IF  C$  = 

•G"  OR  C$  =  "g"  THEN  GOTO  390 

380 

IF  C$  = 

"C"  OR  C$  =  "c"  THEN  GOTO  120  ELSE  GOTO  360 

390 

CLS 

400 

F  =  (B  * 

B  -  A  *  C):  IF  F  =  0  THEN  GOTO  900 

410 

XI  =  U  + 

(B  *  E  -  C  *  D)  /  F:  Yl  =  V  +  (A  *  E  -  B  * 

D)  / 

F:  GOSUB  1000 

420 

R2  =  Rl: 

B2  =  Bl  *  180  /  PI:  IF  K  =  1  THEN  B2  =  B2 

-  180 

IF  B2 

<  0  THEN 

B2 

=  360 

+  B2 

430 

T  =  SGN(B)  *  PI  /  4:  IF  A  =  C  THEN  GOTO  450 

440 

T  =  .5  * 

ATN(2  *  B  /  (A  -  C)) 

450 

G  =  (C  * 

COS(T)  *  COS(T)  -  2  *  B  *  COS(T)  *  SIN(T) 

f  A  * 

SIN(T) 

*  SIN(T)) 

/ 

(-F): 

G  =  SQR(G) 

460 

H  =  (C  * 

SIN(T)  *  SIN(T)  +  2  *  B  *  COS(T)  *  SIN(T) 

f  A  * 

COS(T) 

*  COS(T)) 

/ 

(-F): 

H  =  SQR(H) 

470 

IF  H  >=  ( 

5  GOTO  490 

480 

Z  =  H:  H 

=G:G=Z:T=T+PI/2 

490  CLS 

500  PRINT  "bearing  estimate  =  ";  B2 

510  PRINT  "range  estimate  =  ";  R2:  PRINT 

520  IF  P$  =  "N"  THEN  GOTO  550 

530  LPRINT  :  LPRINT  :  LPRINT  "bearing  estimate",  B2 

540  LPRINT  "range  estimate",  R2 

550  INPUT  "confidence  region  by  size  or  probability  (s/p)";  C$ 

560  IF  C$  =  "S"  OR  C$  =  "s"  THEN  GOTO  580 

570  IF  C$  =  "P"  OR  C$  =  "p"  THEN  GOTO  650  ELSE  GOTO  550 

580  INPUT  "confidence  region  size";  S 

590  IF  S  <=  0  THEN  GOTO  580 

600  CP  =  1  -  EXP(-S  *  S  /  2) 

610  PRINT  :  PRINT  "containment  probability  =  ";  CP 

620  IF  P$  =  "N"  THEN  GOTO  720 

630  LPRINT  :  LPRINT  "confidence  region  size",  S 

640  LPRINT  "containment  probability",  CP:  GOTO  720 

650  INPUT  "containment  probability  (confidence)";  CP 

660  IF  CP  >=  1  OR  CP  <=  0  THEN  GOTO  650 

670  S  =  SQR(-2  *  LOG(l  -  CP)) 


680  PRINT 
690  IF  P$ 
700  LPRINT 
710  LPRINT 
720  X  =  2 
730  PRINT 
740  N  =  T 
750  PRINT 
760  Y  =  2 
770  PRINT 
780  PRINT 
790  IF  P$ 
800  LPRINT 


PRINT  "confidence  region  size  =  ";  S 

"N"  THEN  GOTO  720 
LPRINT  "containment  probability",  CP 

confidence  region  size",  S 

S  *  H 

length  =  ";  X 

IF  N  <  0  THEN  N  =  N 

direction  = 


major  axis 

*  180  /  PI: 
"major  axis 

*  S  *  G 

"minor  axis  length  = 
"containment  ellipse 
=  "N"  THEN  GOTO  840 
major  axis  length" 


180 


N 


";  Y 
area 


PI 


N 


810  LPRINT  "major  axis  direction", 

820  LPRINT  "minor  axis  length",  Y 

830  LPRINT  "containment  ellipse  area". 


PI 


Y  /  4 


Y  /  4 


840  PRINT 
850  IF  C$ 
860  IF  C$  =  "Q" 
870  CLS  :  INPUT 


PRINT  :  INPUT  "continue  or  quit  (c/q)";  C$ 
"C"  OR  C$  =  "c"  THEN  GOTO  870 

OR  C$  =  "q"  THEN  END  ELSE  GOTO  840 

"continue  data  input  or  generate  an  estimate  (c/g)' 

'g"  THEN  PRINT  :  GOTO  500 
OR  C$  =  "c"  THEN  GOTO  120  ELSE  GOTO  430 
900  PRINT  "no  solution" 
910  IF  P$  =  "N"  THEN  GOTO  930 
920  LPRINT  :  LPRINT  "no  solution" 
930  END 
940  XI  =  U  -  R  *  SIN(Q):  Yl  =  V  -  R 


C$ 


880  IF  C$ 
890  IF  C$ 


=  "G"  OR  C$ 
=  "C" 


950  W  =  P  -  Bl:  L 


Rl 


COS(Q):  GOSUB  1000 
BE:  IF  L  =  0  THEN  GOTO  900 
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960  6  =  COS(Bl)  /  L:  H  =  SIN(Bl)  /  L 

970  IF  W  >=  PI  THEN  W  =  W  -  2  *  PI:  GOTO  990 

980  IF  W  <=  -PI  THEN  W  =  W  +  2  *  PI 

990W=W/BE:A=G*G+A:B=G*H+B:C=H*H+C:D=W*G+D:E=U*H+E:  RETURN 

1000  Rl  =  SQR(X1  *  XI  +  Yl  *  Yl):  IF  Rl  =  0  THEN  Bl  =  0:  RETURN 

1010  IF  ABS(X1  /  Rl)  =  1  THEN  Ml  =  SGN(Xl)  *  PI  /  2  ELSE  Ml  =  ATN(X1  /  Rl  /  SQR ( 1  -  XI  *  XI  /  Rl  /  Rl)) 

1020  IF  ABS(Y1  /  Rl)  =  1  THEN  Bl  =  (1  -  SGN(Yl))  *  PI  /  2  ELSE  Bl  =  PI  /  2  -  ATN(Y1  /  Rl  /  SQR(1  -  Yl  *  Yl  / 

Rl  /  Rl)) 

1030  IF  Ml  <  0  THEN  Bl  =  2  *  PI  -  Bl 

1040  RETURN 
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